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ABSTRACT 


We characterize the radial density, metallicity and flattening profile of the 
Milky Way’s stellar halo, based on the large sample of spectroscopically con¬ 
firmed giant stars from SDSS/SEGUE-2 (Xue et al. 2014), spanning galactocen- 
tric radii 10 kpc < roc < 80 kpc. After excising stars that were algorithmically 
attributed to apparent halo substructure (including the Sagittarius stream) the 
sample has 1757 K giants, with a typical metallicity precision of 0.2 dex and a 
mean distance accuracy of 16%. Compared to blue horizontal branch stars or 
RR Lyrae variables, giants are more readily understood tracers of the overall 
halo star population, with less bias in age or metallicity. The well-characterized 
selection function of the sample enables forward modeling of those data, based 
on ellipsoidal stellar density models, z), with Einasto profiles and (broken) 

power laws for their radial dependence, combined with a model for the metallicity 
gradient and the flattening profile. Among models with constant flattening, these 
data are reasonably well fit by an Einasto profile of n — 3.1 ±0.5 with an effective 
radius r e g = 15 ± 2 kpc and a flattening of q = 0.7 ± 0.02; or comparably well 
by an equally flattened broken power law, with radial slopes of a vn = 2.1 ± 0.3 
and a out = 3.8 ± 0.1, with a break radius of ri ireo j. = 18 ± 1 kpc; this is largely 
consistent with earlier work. We find a modest but significant metallicity gra¬ 
dient within the ‘outer’ stellar halo, [Fe/H] decreasing outward. If we allow for 
a variable flattening q = f(i'Gc), we find the distribution of halo giants to be 
considerably more flattened at small radii, g(10 kpc) = 0.55 ± 0.02, compared 
to q(> 30kpc) = 0.8 ± 0.03. Remarkably, the data are then very well fit by a 
single power law with index of 4.2 ± 0.1 on the variable r q = R 2 + (z/q(r)) 2 . 
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In this simple and better-fitting model, there is a break in flattening at ~ 20 kpc, 
instead of a break in the radial density function. While different parameteriza- 
tions of the radial profile vary in their parameters, their implied density gradient, 
<91n v*/dlnr, is stable along a direction intermediate between major and minor 
axis; this gradient is crucial in any dynamical modeling that uses halo stars as 
tracers. 

Subject headings: galaxies: individual (Milky Way) - Galaxy: halo - Galaxy: 
stellar content - stars: K giants 


1. Introduction 


The Milky Way’s extended stellar halo contains only a small fraction (< 1%) of the 
Galaxy’s stars but is an important diagnostic of its formation and dark matter distribution. 
The position-kinematics-abundance substructure in the stellar halo reflects the Galaxy’s 
formation history, whether halo stars were born in situ or are disrupted satellite debris. By 
now, individual stars are the by far the largest sample of kinematic tracers, with 10 kpc < 
vqq < 100 kpc (as opposed to globular clusters or satellite galaxies), and are hence the best 
tracers to determine the mass profile of Milky Way’s dark matter halo in this radial range. 
It is obvious that good kinematic tracer samples should be large and cover a wide radial 
range with accurate individual distances. However, beyond this, the spatial distribution 
of the tracers — in particular their radial profile — must be well understood to use such 
tracers for dynamical inferences. This point is perhaps clearest when considering the |Jeans 


(1915) equation, even in its simplest, spherical and isotropic version: the tracer density 


profile, z/*(r), in particular its logarithmic radial derivative, dlniq/cfinr, affects the inferred 
enclosed mass, M(< r), almost linearly. If we do not know the local power law exponent to 
better than, say, 25%, we cannot infer the mass to better than 25% irrespective of the size 
and quality of the kinematic sample. To a somewhat lesser extent, the inferred mass also 
depends on the flattening of the tracer population. Yet, at present there is little consensus 
on the shape and the radial profile of the stellar halo. 

The most straightforward way to quantify the stellar halo distribution is via star counts. 
Because this method requires large samples of well-understood completeness, it is often 
applied to photometric catalogs. Early studies adopted star counts to analyze globular 
clusters (Harris||1976), RR Lyrae variables (|Hawkins 1984; Wetterer fe McGraw||1996; Vivas 


& Zinn 2006), blue horizontal branch (BHB) stars (Sommer-Larsen 1987), a combination 
of BHBs and RR Lyraes (Preston et al. 1991), a star sample near the north galactic pole 


(Soubiran 1993), or K dwarfs (Gould et al. 1998) and found that the stellar halo is well 
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fitted by a single power law (SPL, v ~ (r Gc)~ a ) with index a = 3 — 3.5 and flattening of 
q = 0.5 — 1. However, Saha (1985]) found that RR Lyrae are well described by a broken 
power law (BPL) with a ~ 3 out to 25 kpc, and a ~ 5 beyond 25 kpc. 

These earlier studies were based on a few hundred objects at most. In recent years, with 
the development of large sky surveys, the sizes of the photometric samples have expanded 


by a factor of more than 10. Robin et al. (2000) used a wide set of deep star counts in a 


pencil-beam survey at high and intermediate Galactic latitudes to model the density profile 
and found the best-fit density profile with a flattening of 0.76 and a power index of 2.44. 


Siegel et al. (2002) found that 70,000 stars in seven Kapteyn selected areas are consistent 


with a power law density with index of 2.75 and flattening of 0.6. However, both Robin 


et al. (2000) and Siegel et al. (2002) used galaxy models to work out the contamination of 


star counts by the disk population. This technique is not ideal for the lowest-density stellar 
component in the Galaxy because any errors in the disk or thick disk will overwhelm the halo 


results. Morrison et al. (2000) used halo stars near the main-sequence (MS) turnoff from the 
“SPAGHETTI” survey to map the Galactic halo and found that the halo density law over 
Galactocentric radii of 5-20 kpc and z -heights of 2-15 kpc followed a flattened power law 
halo with the flattening of 0.6 and a power index of 3. Bell et al. (2008) used ~ 4 million 


color-selected MS turnoff stars from the fifth data release of the Sloan Digital Sky Survey 
(SDSS) up to 40 kpc to find a “best-fit” oblateness of the stellar halo of 0.5 — 0.8, and the 
density profile of the stellar halo is approximately described by a power law with index of 
2 — 4. Subsequently, Sesar et al. (2011) used 27, 544 near-turnoff MS stars out to ~ 35 kpc 


selected from the Canada-France-Hawaii Telescope Legacy Survey to find a flattening of 
the stellar halo of 0.7 and the density distribution to be consistent with a BPL with an 
inner slope of 2.62 and an outer slope of 3.8 at the break radius of 28 kpc, or an equally 


good Einasto profile (Einasto 1965) with a concentration index of 2.2 and effective radius of 


22.2 kpc. 

Deason et al. 

(2011 

analyzed 

~ 20, 000 A-type photometric stars selected from 

the SDSS data release 8 

(Ahn et al. 

2012 

and obtained a best-fitting BPL density with an 


inner slope of 2.3 and an outer slope of 4.6, with a break radius at 27 kpc and a constant 


flattening of 0.6. Subsequently, Deason et al. (2014) found a very steep outer halo profile 
with a power law of r -6 beyond 50 kpc, and yet steeper slopes of a = 6 — 10 at larger radii. 

In addition, several pieces of work point to variations of the stellar halo flattening with 


radius. Preston et al. (1991) found that the density distribution of RR Lyrae follows a power 


law with a ~ 3.2, together with a variable flattening changing linearly from 0.54 at center 


to 1 at 20kpc. Subsequent work (Deason et al. 2011; Sesar et al. 2011) found evidence for 


flattening, but no evidence for a change with radius. 


Spectroscopic maps of the stellar halo beyond ~ 20 kpc in practice require luminous 
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post-MS stars, as turnoff or other MS stars are too faint for current wide-held spectrographs. 
RR Lyrae and BHB stars have repeatedly been used as tracers to study the halo density 
profile, because they have precise distances and are bright enough to be observed at radii out 
to ~ 100 kpc (Xue et al. 2008, 2011; Sesar et al. 2010; Deason et al. 2011, 2014). Yet, such 


stars, most prevalent in particularly old and metal-poor populations (Bell et al. 2010; Xue 


et al. 2011), are known to have a different structure and profile from the metal-poor red giant 


branch (RGB) stars (K-giants). Moreover, Preston et al. (1991) and Morrison et al. (2009) 
claimed that BHB stars in the inner halo had a different radial profile from RR Lyraes. To 
this end, it is crucial to construct the halo shape and radial profile of the stellar halo in K 
giants. Xue et al. (2014) presented a catalog of K giants with unbiased distance estimates 
good to 16% , metallicities, velocities, and photometric information, drawn from the Sloan 
Extension for Galactic Understanding and Exploration (SEGUE; Yanny et al. 2009), which 
contains ~ 150 stars beyond 60 kpc. SEGUE was focused on the halo in its targeting and 
has two phases, SEGUE-1 and SEGUE-2. The target selection of K giants changed a lot 
during SEGUE-1, but SEGUE-2 adopted a unified K-giant selection, so one can understand 
and model their selection function well for K giants observed in SEGUE-2. Owing to the 
well-understood selection function, it is possible to determine the halo profile and shape of 
these tracers, which is the main goal of the present paper. Specifically, we set out to describe 
the stellar halo distribution, presuming that the density is stratified on (oblate) spheroids, 
with a radial profile from 10 to 80 kpc that can be characterized by simple functional forms 
(either an Einasto profile or BPL). We also explore the metallicity dependence of the shape 
and radial profile of the stellar halo. 


In the next section, we lay out the properties and the selection function of the SEGUE K 
giants. In §3, we present the method of fitting a series of parameterized models to SEGUE- 
2 K giants, explicitly and rigorously considering the selection function. This step is key 
in obtaining accurate radial profiles and metallicity gradient model in the meantime. The 
results for the stellar halo’s radial profile and flattening are presented in §4, along with the 
metallicity gradient in the stellar halo. Finally, §5 discusses the comparison between our 
results and previous work, as well as implications for dynamical models. 


2. SEGUE-2 K giants and their selection function 


The SDSS (SDSS; York et al. 2000) is an imaging and spectroscopic survey covering 


roughly a quarter of the sky, which has both ugriz imaging ( 

Fukugita et al. 

1996; 

Gunn 

et al. 1998; Stoughton et al. 2002; 

Pier et al. 2003; 

Eisenstein et al. 2011 

and low-resolution 


spectra (A/AA ~ 2000). SEGUE is one of the key projects and has two phases, SEGUE- 






















































- 5 - 


1 and SEGUE-2, which aim to explore the nature of stellar populations from 0.5 to 100 


kpc (Yanny et al. 2009, and Rockosi et al. in prep.). SEGUE-2 spectroscopically observed 
approximately 120,000 stars, focusing on the stellar halo of the Galaxy. 

To understand the underlying spatial distribution of the K giants on the basis of this 
spatially incomplete sample, we need to understand (and account for) the probability that a 
star of a given luminosity, color, and metallicity ends up in the sample, given its direction and 
distance. Spectroscopic surveys of the Milky Way are inevitably affected by such selection 


effects (see Rix & Bovy 2013), often referred to as “selection biases.” They arise from a set 


of objective and repeatable decisions of what to observe, necessitated by the survey design. 
In particular, only a small fraction of the sky was covered by SEGUE plates, and for most 
plates only a fraction of stars that satisfy the photometric selection criteria could be targeted 
with fibers. Finally, not all targeted stars yield spectra good enough to result in a catalog 
entry, i.e. had signal-to-noise ratio (S/N) high enough to verify that they are giants and 


yield a metallicity measurement. Bovy et al. (2012) and Rix & Bovy (2013) spelled out how 


to incorporate this selection function in fitting a parameterized model for the stellar density, 
and we follow their approach in this and the next Section. 

Both the SEGUE-1 and SEGUE-2 surveys targeted halo giant star candidates, using 
a variety of photometric and proper-motion cuts. About 90% of the final K-giant sample 
came from objects observed as / -color K-giant targets. The /-color (Lenz et al. 1998) is a 
photometric metallicity indicator for stars in the color range 0.5 < (g — r) 0 < 0.8, designed 
to select metal-poor K giant^j As mentioned in fJT] only the selection function for K giants 
observed in SEGUE-2 can be understood well because SEGUE-2 adopted a consistent color- 
magnitude cut to select K giants throughout its entire survey: 15.5 < g 0 < 18.5, r 0 > 15, 
0.7 < (u — g) o < 3, 0.5 < (g — r) 0 < 0.8, 0.1 < (r — *) 0 < 0.6, and l —color > 0.09. Therefore, 
we restrict our analysis to this category. We also require that /-color K giant candidates have 
good proper-motion measurements <11 mas yr _1 . Since broadband photometry is a poor 
MS versus giant discriminator, not all stars targeted under the above criteria will be giants. 
The subsequent identification of K giants is based solely on their spectroscopic properties. 


As described in Xue et al. (2014), this requires spectra that have good Mg index and stellar 


atmospheric parameters determined by the SEGUE Stellar Parameter Pipeline (Lee et al. 
2008ajb, 2011) but have no strong G band. 


To understand the selection function, we must compare the color-magnitude distribution 
of these spectroscopically confirmed K giants with the analogous distribution of all possible 
photometric / -color K giant candidates. Figure [I] shows these two distributions (as con- 


1 https://www.sdss3.org/dr9/algorithms/segue_target_selection.php^S2_table. 



























6 


tours and grayscale, respectively), summed over all SEGUE-2 plates. As the marginalized 
histograms at the sides of the panels show, these two distributions are nearly indistinguish¬ 
able: the chance of a photometric candidate being confirmed as a K giant is independent of 
color and magnitude. This simplifies the subsequent analysis and is testament to SEGUE-2’s 
consistency of target selection, targeting, and spectral analysis. 

However, while the selection function is constant with apparent magnitudes and dered- 
dened colors, it varies from plate to plate, in particular with the Galactic latitude of the 
plate. Given the pencil-beam nature of the SEGUE survey, it makes sense to specify the 
selection fraction plate by plate. For each plate we define the number of spectroscopically 
confirmed K giant as N spec , and the number of /-color K-giant candidates in the plate (both 
those that were targeted and those that were not) as N p h ot . Thus, the plate-dependent 
selection function (shown as Figure [ 2 ]) is given by 

S (plate) = (1) 

™ phot 


As we want to analyze the spatial distribution, we also restrict our sample to SEGUE-2 
/-color K giants that have reliable distance estimates from Xue et al. (2014). To eliminate 
the contamination from the disk component, we cull K giants with [Fe/H] > —1.2 and 
|z| < 4 kpc, which leads to a final sample of 2413 1-color K giants. Figure [3] illustrates the 
basic sample properties: its sky coverage and spatial distribution (without accounting for 


the selection function), and the distribution of metallicity against distance from Xue et al. 


(2014). The stars’ distribution reflects the pencil-beam pattern of the SEGUE survey; their 
galactocentric distances range from 7 to 85 kpc; their mean metallicity is —1.75 dex, with 
some being as metal-poor as —3.5. The reader should note that these metallicities are not 
typical of the halo as a whole because of our choice to exclude all stars with [Fe/H] > —1.2. 


3. Modeling the Stellar Density and Metallicity Distribution in the Halo 


In this section we lay out a forward-modeling approach to describe the spatial and 
metallicity distribution of the stellar halo with a set of flexible but ultimately smooth and 
symmetric functions. Both the modeling practicalities and the astrophysics require that we 
fit the spatial distribution and the abundance distribution simultaneously. We defer to the 
§4 the question of how sensitively such a “smooth” description depends on the question of 
including or excluding stars that are presumably members of recognizable substructure (see 


Belokurov et al. 

2006; 

Bell et al. 

2008) 
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3.1. A Parameterized Model for the Observables 


We presume that the stellar halo distribution can be sensibly approximated by a spheroidal 
distribution with a parameterized radial profile, allowing for radial variations in the metallic- 
ity distribution, ^*(x, y, z, [Fe/H]) = p([Fe/H] | r q , pF e ) x 14 (r q | pn) . The combination of pn 
and ppe denotes the halo parameters , and r q = \[W‘ + (z/g(r)) 2 is the basic galactocentric 
radial coordinate, given the flattening q(r). 


In this section, we spell out a straightforward and rigorous approach to determine the 
posterior probability distributions for halo parameters in light of the above data, our knowl¬ 
edge of the SEGUE-2 selection function, and well-established astrophysical priors on the 
luminosity function of giant stars. The number of halo parameters depends on the com¬ 
plexity of the model stellar halo distribution. This approach essentially follows [Bovy et al. 


(2012 

and 

Rix & Bovy 

(2013) 


Since we already have good estimates of the distance modulus VA4 and the r-band abso¬ 
lute magnitude M r for all objects in the sample (Xue et al. 2014), we treat (V-,, M r , [Fe/H]) 
as the observables defining V Q = ("DAT l, b ), rather than (rn, c, [Fe/H], l, b), as it makes the 
fitting formalism more intuitive. However, the r-band apparent magnitude m and (dered- 
dened) color c (referring in particular to g — r) will appear explicitly in the observational 
selection function. We denote the angular selection function as S(l,b) = ,S'(plate(l. bj) (see 
Eq.(l)), the magnitude-color selection function by S(m(V A4, M r ), c(M r , [Fe/H])), the metal- 
licity selection function by S([Fe/H]), additional spatial cuts by S(V Q ), and the luminosity 
selection function as S ( M r | [Fe/H]), expressed in terms of the “observables” above. We 
denote the prior external information on the distribution of absolute magnitude on the gi¬ 
ant branch as p(M r ); for old and metal-poor populations this is well established through 
cluster luminosity functions (see, e.g., Xue et al. 2014). The luminosity selection function 
S(M r | [Fe/H]) will appear explicitly, because stars low on the giant branch were removed 
(depending on [Fe/H] ) in Xue et al. (2014) to avoid confusion between RGB and red clump 
(RC) stars. Figure [4] shows the limits of M r for a given [Fe/H], In short, we need to spell out 
and then quantify all terms that matter for predicting the rate function , i.e. the expected 
frequency or probability of finding a sample member with a given VAi, l, b, M r , [Fe/H] if a 
halo with pu and py e were true. 


Given p H and py,. , the expected rate function for finding a star with ('D-., M r , [Fe/H]) is 


then 
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A(P Q ,M r , [Fe/H] | p H ,PFe) = 

|J(x,y,z;D 0 )| x i/*(X> 0 | p H ) x p([Fe/H] | £> 0 ,p Fe ) x p(M r )x 
S{V 0 ) x S(l,b) x S([Fe/H]) x S(M r | [Fe/H]) x s(m(VM, M r ),c(M r , [Fe/H])). 

( 2 ) 


This rate function is simply a concise way of specifying what set of observations we 
expect, given an assumed halo model, including all observational selection effects and perti¬ 
nent prior information. The price for putting this in one place is that the rate function has 
quite a number of distinct terms. In the next two subsections, we spell out and discuss these 
terms. To make the rate function a probability, it must be normalized for every new set of 
trial model parameters (pu , p Fe ); this is the only time-consuming step in such modeling. 


3.2. Models for the Spatial and [Fe/H] Distribution of the Stellar Halo 

Following a number of previous studies, we presume that the overall radial density profile 
of the halo can be described by an Einasto profile or by (multiply) broken power law, with 
the density stratified on surfaces of constant r q in all cases; in addition, we devise a simple 
model for radial [Fe/H] variations. 


3.2.1. Einasto profile 


The Einasto profile (Einasto 1965) is the 3D analog to the Sersic profile (Sersic 1963) 


for surface brightnesses and has been used to describe the halo density distribution (Merritt 


et al. 

2006; 

Deason et al. 2011; 

Sesar et al. 

2011 

): 

v*(r q ) = vo exp j -d n (r ? /r eff ) 1/n - 1 

}• 

(3) 


where v 0 is the (here irrelevant) normalization, r e g is the effective radius , n is the concen¬ 
tration index, and d n ~ 3n — 1/3 + 0.0079/n, for n > 0.5. The free parameters of an Einasto 
profile with a constant flattening are pn = (r e ff,n, q). 
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3.2.2. Broken Power-law Profiles 


Broken (or even multiply broken) power laws (BPLs) are another family of functional 
forms that has been used extensively to described the radial profile of the Galactic stellar 
halo (Saha 1985; Deason et al. 2011, 2014; Sesar et al. 2011). In most cases the change in 
the power law index, dlnz^/dlnr, has been taken to be a step function. We adopt 


*4 ( r q) 


u 0 xr 


U) r q “ i: 

(ctout Gfiri) 

break 


X r 


C^out 

q ’ 


I*q ^ I*break 


(4) 


In addition to the flattening and the (irrelevant) normalization, a (singly) broken power law 
has three parameters: a lTl , a ouf , r^eak . Of course, this profile family encompasses an SPL, 
and it can be generalized to a multiple-broken power law (twice-broken power law, TPL) by 
introducing an additional pair of («, rbreak ) ( see 


Deason et al. 


2014 


3.2.3. Halo Flattening 


While Preston et al. (1991) found evidence for a decrease of flattening with increasing 


radius, others did not (Deason et al. 2011; Sesar et al. 2011). As there is evidence that at 


least the innermost part of the halo is quite flattened (Carollo et al. 2010), we explore how 


our sample can inform us about radial variations in the flattening of the stellar halo beyond 
rcc =10 kpc. To date no particular functional form to parameterize the possible variation 
of halo flattening has been established; therefore, we consider the functional form for q(r) 


as: 


9 (r) = qoo - (hoc - qo) exp 1 - 


Vr 2 + : 

ro 


(5) 


where qo is the flattening at center, changing to q^ at large radii, with r 0 is the (here 
exponential) scale radius, over which the change of flattening occurs. 


3.2.4- Radial Variations in the Metallicity Distribution 


The existence of a halo metallicity gradient is currently controversial (see, e.g., Carollo 


et al. 2007; Schonrich et al. 2011; Fernandez-Alvar et al. 2015). This is partially due to the 


lack of in situ tracers for which accurate abundance measures are possible with low-resolution 
spectroscopy, and partially due to the lack of attention to the effects of the complex SEGUE 
selection function. This is the first paper that has applied a forward-modeling approach to 
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the K giants in the halo, addressing both of these concerns. We choose a mixture model for 
the halo metallicity distribution as follows. The metal-rich component can be described by 
a Gaussian with a mean at —1.4 dex and a sigma of 0.2 dex, and the metal-poor component 
follows a Gaussian with a mean at — 2.1 dex and a sigma of 0.35 dex. Therefore, we suppose 
the metallicity distribution , p([Fe/H] | V Q ,p Fe ) as a radially varying combination of these two 
metallicity components. The metallicity distribution model is also expressed as a function 
of r q . We choose 

P([Fe/H] |r q , p Fe ) = f 0 x (r q / 20 kpc)^ x £(-1.4,0.2) + (l-f 0 x (r q / 20 kpc)^) x £(- 2 . 1 ,0.35) ( 6 ) 

where r q is the same as defined in Eq. [3j £(mean, dispersion) are Gaussian functions. Besides 
the common parameter q, the metallicity gradient has two other parameters, f 0 and 7 (we 
mark them as pp e ). 


3.3. Selection Effects and Prior Information 


The remaining terms in Eq. 2 incorporate the various selection effects and pieces of 
prior information into the prediction of the rate function and are given in the following. 


The Jacobian term is given by 


\J(x,y,z;T>M)\piate = Opiate ■ In 10/5 • ( 10 ™ 2 kpc) J , 


(7) 


where f2 p i ate = 7 deg 2 . The prior on the absolute magnitude distribution along the giant 


branch (Xue et al. 2014) is given by 

10 o-32 Mr , if M r minj0bs < M r < M r maXj0bs 

0 , otherwise, 


p(M r ) oc 


( 8 ) 


The sample’s metallicity cuts, aimed at eliminating bulge and thick-disk contributions, 
as well as any spurious metallicity determinations, —3.5 < [Fe/H] < — 1 . 2 , are reflected in 


S([Fe/H]) oc 


1. « [Fe/HU,* < [Fe/H] < [Fe/H] maxobs 

0 , otherwise 


( 9 ) 


The spatial cuts to geometrically excise any bulge and thick-disk stars are 


S(D®)) (X 


1, if |z(D 0 )| > 4 kpc 
0 , otherwise 
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( 10 ) 


The reliable distance determination requires minimizing the contamination of RGB stars 
by RC stars, which translates into [Fe/H]-dependent absolute magnitude cuts (Figure |4j) : 


S(M r | [Fe/H]) 


1, if M min ([Fe/H]) < M r < M max ([Fe/H]) for [Fe/H] G [-3.5, -1.2] 

0, otherwise 

( 11 ) 


The set of observed SEGUE plates leads to a rather sparse angular selection function 
in (/,&), illustrated in Figures [2] and |3j 


S (lpla.te! bplate ) 


if in plate 

phot 

0, otherwise, 


( 12 ) 


Finally, the apparent magnitude and color cuts (see Figure [I]) for the SEGUE targeting 
are reflected in 


S(m r (lTM,M r ),c(M r ,[Fe/H])) oc 


1, if m m i n O m r O m max and c m j n <C c O c 
0, otherwise 


(13) 


Taking these terms together and following Bovy et al. (2012), we can now directly 


calculate the likelihood of the data given (pn, pFe) and the rate function: 

Nkg 

^(datajlp^PFe) = c“ Nkg JjA(M ri ,PM i ,[Fe/H] i ,l i ,b i |pH,PFe), 

i=l 


(14) 


where i is the index of each K giant and we use that the data points are identically 
and independently drawn. The normalization c\ is the integral over the volume in the 
(T>A4, 1, b, M r , [Fe/H]) space, 


Nplate f f f 

C A = J J j [Fe/H],l plate ,b plate |p H ,PFe)dM r dDA4d[Fe/H] (15) 

Here we have already performed the dldb integral for each plate. This normalization integral 
is the computationally most expensive part of the parameter estimates. But it can be com¬ 
puted efficiently using Gaussian quadratures, where we adopt 48 x 20 x 20 transformation 
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points in VA4. M r and [Fe/H] space, and where the parameter-independent parts such as the 
Jacobian term, the luminosity prior, and selection functions can be pre-computed on a dense 
grid. We assume the priors on all parameters (pH,PFe) to be flat over the pertinent range; 
therefore, the posterior distribution function ( pdf) of the parameters, p(pH,4>Fe|{data}), is 
proportionate to the likelihood (Eq. [l4j), differing only in its units (‘1/parameters’ vs. 
T/data’). We then vary the (pH - pFe) to sample the parameter pdf using emcee, which im¬ 
plements an efficient Markov chain Monte Carlo techniqu^] (Foreman-Mackey et al. 2013). 


4. Results 


We now present the results of applying the modeling from §3 to the sample of §2. 
We illustrate these results in two ways: first, by showing the joint pdfs of the halo model 
parameters; second, we show what the radial and flattening profiles actually look like, by 
drawing samples from these parameter pdfs. We start out with the simplest model (the 
one with the fewest parameters), an Einasto profile of constant flattening, to illustrate the 
results, present the basic gradient in the [Fe/H] distribution, and to explore the impact of 
excising stars that are manifestly in substructures. We then proceed to include variable 
flattening, and the BPLs as radial profiles. The results are summarized in Table 1. 


Figure [5] illustrates the result of the simplest model we fit to the entire data set, by 
sampling the parameter pdf using Eq. 14 This model has three parameters for zp (r q \ph), 
Ph = {r e ff, n, (]}., and two parameters for the metallicity distribution p([Fe/H] | r q , pp e ) , 
namely pp e = {/o, 7 }. The projections in Figure [5] of the pdf involving only pu are shown in 
gray, while those also involving pp e are shown in blue. This Figure illustrated that the sample 
size and data quality of the sample are sufficient to provide formally very well constrained 
parameters. 


4.1. Radial Gradients in the Halo’s Metallicity Distribution 

Figure [ 5 ] shows that the mix of the intermediate (([Fe/H]} = —1.4) and metal-poor 
(([Fe/H]) = — 2 . 1 ) [Fe/H] components changes with radius (i.e. 7 7 ^ 0): the relative im¬ 
portance of the metal poor component increases toward large radii, by a factor of 1.4 ± 0.1 
from 10 to 65 kpc. Our best-fit metallicity distribution predicts a mean metallicity of -1.8 in 
this radial range. The Figure also shows that the covariances between the fits to the stellar 


2 The user guide for emcee can be found at http://dan.iel.fm/emcee/current/ 
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density and to the [Fe/H] distribution are weak. The resulting [Fe/H] distribution, projected 
into the space of the rate function, is shown in the inset of Figure [5| this inset illustrates 
that the distribution of the actual sample members in the rcc — [Fe/H] plane is plausible. 
Note that nearly all sample stars within rcc = 10 kpc, much of the “inner halo” of Carollo 
, have been eliminated from the fit (see Figure [3]); therefore, our result is the 
detection of an outward metallicity gradient within the “outer” halo. 


et al. (2007 


4.2. The Impact of Halo Substructures on Fitting Smooth Models 


As discussed in §1 , both cosmological models and observations imply that a good 
portion of halo stars, at least beyond 20 kpc, are in substructures. Especially the prominent 
ones, such as the Sagittarius stream and the Virgo overdensity, can and will affect the fits 
of smooth models, as pointed out by Deason et al. (2011). Recently, Janesh et al. (2015) 
used a position-velocity clustering estimator (the 4-distance) in combination with a friends- 
of-friends (FoF) algorithm to identify the stars in the Xue et al. (2014) K-giant sample 
belonging to substructures; they found that 27% of the K giants are clearly associated with 
the Sagittarius streams, the Orphan streams, the Cetus Polar stream, and other unknown 
substructures. The results of Janesh et al. (2015) provide a straightforward algorithmic way 
of excising much of the manifest substructure from the sample. For the Sagittarius stream 
(Belokurov et al. 2014) we know the expected position-velocity distribution quite well, which 
suggests that six of the most distant sample members are likely members of Sagittarius’s 
trailing arm. These he at rcc >60 kpc, 120° < L e < 140°, |B 0 | < 13° and cluster tightly 
around Vg iSr = 100 kms -1 , which also eliminates them from the original sample of 2413 
l — color K giants, leaving 1757 l — color K giants. 


We then repeat the fitting of the same model as illustrated in Figure [5| and find - un¬ 
surprisingly - significant differences (Figure [ 6 ]): an Einasto profile that is more concentrated, 
has a smaller effective radius and is more flattened, n = 3.1 ± 0.5, r e g = 15 ± 2kpc, and 
q = 0.7 ± 0 . 02 . The metallicity gradient, however, remains very similar with fo = 0.6 ± 0.01 
and 7 = — 0.2 ±0.04. Such a model already provides quite a good fit to the data, as Figure [7] 
shows: the model prediction for the VAi distribution of the sample members, averaged over 
all directions and metallicities and accounting for all selection effects, matches the actual 
VAi distribution quite well. This Figure corresponds to the top-right inset in Figure | 6 j just 
marginalized over metallicities. 
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4.3. Different Functional Forms for the Radial Profile 


Following previous approaches, we next explore different functional forms for the radial 
profile, assuming radially constant flattening. 

Figure [8] shows the results analogous to Fig. [6j but for a BPL: at a best-fit, constant 
flattening of q = 0.7 ± 0.02 we find an inner slope of = 2.1 ± 0.2, distinctly different 
from the outer slope of a ou t = 3.8 ±0.1, with a break radius of rb rea k = 18 ± 1 kpc. The 
flattening q is consistent with that of the best-fit Einasto profile. Again, also with BPL 
profile p(VAi, [Fe/H]) remains similar and is consistent with the observations. Bayesian 


information criterion (BIC; Schwarz 1978, ;also called the Schwarz criterion) is a criterion 
for model selection. It is defined as — 2 In 2zf max ± N p In N ( ] ata , where 2zf max is the maximized 
likelihood of the model, N p is the number of parameters in the model, and A^ata is the 
sample size. BIC takes into account both the statistical goodness of fit and the number of 
parameters that have to be estimated to achieve this particular degree of fit, by imposing a 
penalty for increasing the number of parameters. The model with the lowest BIC is preferred. 
As the values of ABIC in Table 1 indicate, the Einasto and BPL radial profiles fit the data 
comparably well. A twice-broken power law (TPL), with the break radii held fixed at the 
radii suggested by Deason et al. (2014), leads to a comparably good fit (see Fig. [9] and 
Table 1). Compared to Deason et al. (2014), the fit is consistent within 65 kpc; it shows no 
evidence for a steep drop beyond, but our present sample is vary sparse at these distances. 
For reference, we also fit the data with an SPL of constant flattening, but on the basis of its 
ABIC it can be ruled out compared to the BPL and TPL models. Our results confirm that 
if constant flattening is assumed, the halo profile must be described by a radial profile break 
at ~ 20 kpc. 


4.4. Radial variation of halo flattening 


We now proceed to allow more model flexibility by allowing the flattening q{r ) to vary 
with radius, according to Eq. [5| we fit both the Einasto profile and power laws (BPL and 
SPL) to the data. Allowing for flattening variations makes for distinctly better fits to the 
data, as quantified by A BIC (Table 1). These fits imply a strong variation of the halo 


flattening with radius, as the pdfs for an Einasto profile in Fig. [10] illustrate: while at large 
radii q^ ~ 0.8, the flattening at asymptotically small radii becomes formally very strong, 
q 0 ~ 0.2. Given that the characteristic radius at which this flattening change occurs is 
r qo ~ 6 kpc, and hence smaller than the minimal radius of all data points at 10 kpc, the 
actual flattening changed across the radial range constrained by the data is more modest, 
as Fig. [T2] illustrates: q changes from 0.55 ± 0.02 at 10 kpc to 0.8 ± 0.03 at large radii. 
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Allowing for variable flattening also changes the Einasto parameters quite drastically: the 
formal effective radius becomes very small, and there is a strong covariance between r e g 
and concentration n. However, the actual predictions for the slope of the radial profile, 
<91n u^/dlnr, in the radial regime constrained by the data are quite similar between models 

This just illustrates to 


with constant and variable flattening, as Figs. [13] and 14 show, 
compare model fits drawn from different functional form families in their projection in the 
space of observables, not just in the space of the parameters. The analogous pdjs of the 
analogous fit for a BPL look at first sight bewildering. However, it reveals that the implied 
flattening profile is the same as for an Einasto fit. The unconstrained pdjs on the break radius 
simply reflect the fact that the inner and outer power law slopes become indistinguishable, 
which is equivalent to an SPL. Therefore, we turn to fit the data to an SPL with a varying 
flattening (see Figure 0 Remarkably, we find a very good fit with an SPL in r q when 
allowing for variable flattening, because this model has the minimum BIC, as illustrated by 
ABIC = 0 in Table 1. The flattening profile is the same as for an Einasto fit, as also shown 


in Figure [12] there is effectively a break in the flattening profile at 15-20 kpc. Figure 13 


shows that the actual local density slope along an intermediate axis for this SPL is similar 
to the other fits. Taken together, this implies in the context of these density models that 
the stellar halo density is best and simplest described (10-80 kpc) by an SPL in the variable 
r q , with the break in the radial profile at 20 kpc replaced by a break in the flattening at a 
comparable radius (see Figure [l2]). 


5. Discussion and Summary 

Using a large spectroscopically confirmed sample of K giants, well-understood popu¬ 
lation tracers of the stellar halo, we have attempted to characterize the radial profile, the 
flattening profile, and the metallicity profile of the ‘outer halo’ from 10 to 80 kpc; we have 
done so after excising algorithmically identified members of distinct halo substructures from 
the sample. On the one hand, this analysis has reemphasised that the choice of the func¬ 
tional forms for the density models matters in casting the results; on the other hand, we 
could show that three robust aspects emerge: the profile of the local density slope, 
when taken along an axis intermediate between major and minor axis, is consistent among 
all functional profile forms we explored; the stellar halo distribution appears distinctly flatter 
within 20 kpc; and there is a slight, but significant, outward metallicity decrease within in 
the ‘outer halo.’ We also find a break in the properties of the halo density distribution at 
~ 20 kpc. This can be attributed to a break in the radial power law, as in a number of 
previous analyses, but we find that it foremost reflects a break in the halo flattening: the 
data are well fit by an SPL in the spheroidal coordinate r q , with a distinct change in q(r) at 
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a radius of r^j 20 kpc. 


We illustrate the first point about the local slope in Figure [l 3 j As mentioned in the 
Introduction, this slope plays a particular role in dynamical modeling, perhaps most easily 


seen in the case of axisymmetric Jeans equation modeling. Figure 13 compares the radial 
profile functions between the models with constant or variable flattening: all models have 
consistent radial profile slopes, within 65 kpc. There is remarkable qualitative 

agreement, in the sense that these slopes vary consistently with radius, within the constraints 
of the functional form imposed. In all cases does the halo profile steepen between 10 and 
65 kpc. Clearly, the fits with fixed q differ more from one another and from the q(r) fits, 
as they must through their more stringent functional form constraints. Allowing for a yet 
steeper drop beyond 65 kpc with a TPL profile leads to indistinguishable results. The largest 
discrepancies are between power laws and Einasto profiles beyond 65 kpc; we attribute that to 
the Einasto profile being constrained by the abundance of data at r gc < 65 kpc (see bottom 
panel of Fig. [l 3 | , which inevitably leads to at large radii. A look at Table 1 shows, 
however, that all of these models make the data appear comparably likely; the likelihood 
differences are significant, but not drastic (see also Fig. [7j) . In particular, the Einasto and 
BPL profiles lead to near-identical data likelihoods; the model fits with the flattening forced 
to be constant make the data significantly less likely. Remarkably, the SPL with variable 
flattening has nearly the same along the intermediate axis, where the seeming radial 
change in power law index is simply a reflection of the change in the radial coordinate r q , 
attributable to q(r). 


We also compare these findings to other work, in particular that of Deason et al. (2011 


2014). Figure 14 shows that the halo profile, as traced in this analysis via K giants, is 


consistent with the findings of Deason et al. (2011) within 65 kpc. To follow up the findings 


of Deason et al. (2014), a steep drop in the density of BHB stars beyond 65 kpc, we specifically 


fit a TPL with parameters (aq, 012 , ol 3 , q). For comparison, we held the break radii fixed at 
r breaki = 18 kpc and rb r eak 2 = 65 kpc. Yet, as Figure [M] reveals, the halo slopes in BHB stars 
and found here in K giants are formally inconsistent beyond 65 kpc. Yet, there are only 7 
stars beyond 65 kpc, so the paucity of distant K giants precludes a more stringent comparison. 
Besides, another two avenues are conceivable: first, we know that the substructure differs 
between BHB stars and giants and MS turnoff stars (Bell et a~L|2008 ; Xue et al.|2011; Janesh 


et al. 2015). Second and related, it is conceivable than many of the stars at large radii in 


the present K-giant sample could be Sagittarius stream members, even though we made an 
attempt to remove such stars from the sample. 


In summary, we believe that the present study has brought forward a number of new 
aspects: First, we have worked out a forward modeling of the spectroscopic data that has not 
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been applied previously in this context. We believe that, in particular for tracers that are 
not standard candles like BHB stars, such an approach is warranted and powerful. Second, 
we were able to show on this basis (1) that there is an outward metallicity gradient in the 
halo beyond 10 kpc; (2) that the halo is distinctly flatter between 10 and 20 kpc, compared 
to larger radii; this distinct change in flattening suggests that it is more appropriate to think 
of the break in halo profile at 20 kpc as a break in flattening, rather than as a break in 
the radial profile at forced constant flattening; and (3) that there is overall consistency with 
previous analyses when it comes to the dynamically important quantity in the radial 
range 10-65 kpc. 
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Fig. 1.—: Distribution of SEGUE-2 photometric l — color K-giant candidates in color-color 
and color-magnitude space (gray map, black contours, and histogram) and the spectroscopi¬ 
cally sample targeted and successfully analyzed (red contours and histogram). The contours 
contain 68%, 95%, and 99% of the distribution. The spectroscopic sample with parameters 
(including [Fe/H] and VA4 ) is a fair subset of the photometric targets with respect to colors 
and magnitudes. The color and magnitude limits enter the modeling through Eq. 
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Fig. 2.—: SEGUE-2 targeting fraction of photometrically selected of l — color K giants, as 
a function of Galactic coordinates x and y (upper panel), and of Galactocentric coordinates 
x and vertical height z (lower panel). Each line represents a spectroscopic plate, with the 
extent of the line showing the nearest and farthest object targeted. The color indicates the 
fraction of stars that photometrically pass the selection criteria that actually got targeted 
spectroscopically. At high Galactic latitudes, most photometrically eligible targets had spec¬ 
tra taken; at low Galactic latitude, only a modest fraction, which is, however, known for 
each plate and accounted for in the modeling (Eq. 12). 
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r gc (kpc) x (kpc) 


Fig. 3.—: Upper left: sky coverage and the spatial distributions (right panels) of SEGUE-2 
l — color K giants in the sample; they appear as a pencil-beam pattern, due to the nature of 
SEGUE survey. Lower left: distribution of metallicities, along with the galactocentric radii, 
shows that the mean metallicity is about —1.75 dex, and some K giants have metallicities 
of about —3.5. The stars with [Fe/H] > —1.2 and \z\ < 4 kpc are culled because they could 
belong to the disk. The metallicity and distance cuts enter the modeling though Eqs. [9] and 


10; the angular selection function, through Eq. 12 
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Fig. 4.—: Limits of the absolute magnitudes M r changing with the metallicities [Fe/H], 
The bright limits are caused by the metallicity-dependence of the fiducials’ bright tips, while 


the faint limits are due to removal of possible RC stars by Xue et al. (2014) and enter the 
modeling through EqJTTj 
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Fig. 5.—: All one- and two-dimensional projections of the posterior probability distributions 
for the parameters (q, n, r e g, f 0 , 7 ) of the Einasto-profile and metallicity model, shown here 
for the case that no halo substructures are excised. The red lines and squares mark the 
median value of each parameter. The dashed lines show 68 % confidence interval. The black 
sub-figures show the parameter pdf for the spatial density profiles, while the blue sub-figures 
show parameter pdfs for the metallicity distribution model. The top right figure compares the 
data set (black dots) with the model prediction for [Fe/H] — VAi distribution , averaged over 
all directions and accounting for all selection effects. The colored model contours encompass 
99.9%, 95%, 68 %, 50%, and 16% of the normalized model probability, projected into the 
[Fe/H] — PAd-plane. It matches the actual [Fe/H] — VAi distribution well. Of course, 
the model also makes a prediction for the (Z, b ) distribution of stars, but that prediction is 
dominated by the positions of the spectroscopic plates. 
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Fig. 6.—: One- and two-dimensional projections of the posterior probability distributions 
of parameters (q, n, r e g) of the Einasto-profile fit to the sample remaining after excluding 
stars identified by Janesh et al. (2015) as likely members of recognizable halo substructures. 
Analogous to Fig. [5| the basic data-model comparison is shown in the top right panel. The 
comparison to Figure [5] shows that excluding substructures leads to a more concentrated and 
slightly more flattened Einasto profile; the profile’s concentration parameter n is covariant 
with the effective radius parameter, as r e g approaches the inner distance cutoff of the sample 
(10 kpc). 
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Fig. 7.—: Upper panel: comparison between the observed distance distribution and the 
predicted distributions by the best-fitting models, marginalized over all other parameters 
(e.g. [Fe/H], l and b ), and accounting for all selection effects. Lower panel: analogous 

comparison between the observed metallicity distribution and the model predictions (now 
marginalized over DM). All of the best-fitting models fit qualitatively well. Discrepancies 
in the [Fe/H] distribution reflect the fact that we only considered a metallicity distribution 
that could be reflected in radially varying proportions of two Gaussian functions. 
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Fig. 8.—: Same as Figure |6j but for the BPL profile with parameters (q, cq n , a ou t, r break)- 
The red lines and squares mark the median value of each parameter, and the dashed lines 
show 68% confidence interval. This BPL profile fits the data basically as well as the Einasto 
profile (see also Figure [7]). Indeed, the actual prediction for the slope, din /din r , in the 
radial regime constrained by the data is quite similar to that predicted by the Einasto profile 
(Figure [l3j) . 
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Fig. 9.—: Same as Figure]^, but for the TPL profile with parameters (q, oq, aq, 0 : 3 ); the 
break radii are held fixed at 18 and 65 kpc respectively (for comparison with Deason et al. 
(2014J)). The red lines and squares mark the median value of each parameter, and the dashed 
lines show 68 % confidence interval. The best-fitting TPL also fits the data well shown as 
top right panel. The power law index between 18 and 65 kpc, aq, is comparable to the index 
beyond 65 kpc, oq. There is no strong drop beyond 65 kpc. 
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Fig. 10.—: Parameter pdfs for the Einasto profile with varying flattening (see §4.4). The 


fit implies a strong variation of flattening with radius, illustrated in Figure [12} Allowing 
for a varying flattening changes the parameters of the Einasto profile quite strongly: the 
effective radius becomes formally very small, and there is a strong covariance between r e g 
and n. However, the actual prediction for the profile slope, din v*/din r, is quite similar to 


that predicted by the Einasto profile with a constant flattening shown as Figure 13 
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Fig. 11.—: Parameter pdjs for the SPL profile with varying flattening (see §4.4). The red 
lines and squares mark the median value of each parameter, and the dashed lines show the 
68% confidence interval. The flattening profile is the same as for an Einasto fit (see also 
Figure [l2]). The actual predicted din z/*/din r is similar to other models in the radial regime 


constrained by the data shown as Figure 13 
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0 . 9 - 


— exponential q(r) + Einasto 

— exponential q(r) + SPL 
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Fig. 12.—: Radial variation of the stellar halo’s flattening (see £4.4). The coincidence of 
the curves for an Einasto and SPL radial profile illustrates that the flattening profile is 
independent of the radial functional form. Over the observed range, the halo becomes much 
rounder at large radii, from q=0.55 ± 0.02 at 10 kpc to q=0.8 ± 0.03 at large radii. The 
bottom panel shows the actual radial distribution of the sample members used to constrain 
the fit. 
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excl. all substructures 



Fig. 13.—: (Modest) differences between the differing parameterizations of the radial pro¬ 
files: the top five panels show the slope of the stellar density profile dln ^^R^' > (along an 
intermediate axis), for the Einasto profiles with constant (black) and variable (blue) flat¬ 
tening and the BPL profiles with constant (green) flattening, and the SPL with variable 
(magenta) flattening, and the TPL (red) with constant flattening. The lines correspond to 
radial profiles created from 100 samples of the parameter’s pdf. The most likely profile from 
the top panel is repeated below for reference. The two vertical dashed lines indicate the 
two break radii (18 and 65 kpc) of TPL in our fitting. Note that the break radius is in 
spheroidal coordinates r q , while the X-axis here is spherical radius rcc- Despite the fact 
that the Einasto profiles in the two cases (black and blue) have effective radii that differ 
by a factor of two, their radial profiles are very similar within the range constrained by the 
data. Especially, the Einasto profile and BPL with variable flattening (blue and magenta) 
show great consistency. Again, the bottom panel shows the distance distribution of the data. 
Most of the data are within 65 kpc, where all models have very similar predictions for the 
slope, dlni / (R 1 R/y^{2))/dlnR. 
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Fig. 14.—: Comparison of the radial density slope 


dlm/(R,R/x/2) 

dlnR 


found in the present analysis, 


based on K giants, and the slopes found by Deason et al. (2011, 2014) from BHB stars; the 
differently colored lines show our best-fit BPL with constant flattening (green) and best-fit 
TPL with constant flattening (red), best-fit Einasto profiles with constant flattening (black) 


and the best-fit SPL with variable flattening (magenta). The best-fit models of Deason et al. 


(2011, 2014) are shown in black, green, and red but in dashed fines. The two vertical dashed 
lines indicate the two break radii (18 and 65 kpc) of TPL in our fitting. Note that the break 
radius of broken power law is in r q , but here the x-axis is in r gc . The lower panel shows the 
distance distribution of the data. Among the BPLs there are some differences between the 
present fit and the Deason et al. (2011, 2014) results, which used a different population as 


tracers. Notably, the K-giant sample does not point toward a steep halo drop beyond 65 
kpc, although our substructure-cleaned sample contains only a few stars at such large radii. 








































Table 1:: A Summary of Our Best-fitting Models 


Model 

Parameters 

Np 

AlnJ^, ax 

A BIC b 

Excl. all substructures 

[Fe/H]-model 

f 0 = 0.6 ± 0.01, 7 = -0.2 ± 0.05 

2 



+ 





SPL 

a = 3.6 zb 0.1, q = 0.68 ± 0.02 

2 

-35 

54 

Einasto 

n = 3.1 zb 0.5, r e ff = 15 zb 2 kpc 
q = 0.7 zb 0.02 

3 

-16 

25 

BPL 

ain = 2.1 zb 0.3, a Q ut — 3.8 zb 0.1 
r break = 18 zb lkpc, q = 0.7 zb 0.02 

4 

-14 

27 

TPL 

ai = 2.1 zb 0.2, a 2 = 3.8 zb 0.1 
a 3 = 4.8 zb 0.8, q = 0.7 zb 0.02 
r breakl = 18kpc, r break2 = 65kpc 

4 

-13 

26 

[Fe/H]-model 

f 0 = 0.6 ± 0.01, 7 = -0.26 zb 0.06 

2 



+ 





Einasto-q(r) 

n = 7.7 zb 1.5, r e ff = 2 zb 1 kpc 
q 0 = 0.2 zb 0.1, q inf = 0.78 zb 0.05 
r b = 6 zb 1 kpc 

5 

-1 

9 

BPL-q(r) 

a in = 4.2 zb 0.4, a ou t = 4.3 zb 0.3 
r break = 22 zb 23 kpc ,q 0 = 0.2 zb 0.1 
qinf = 0.8 zb 0.03, r b = 6 zb 1 kpc 

6 

0 

14 

SPL-q(r) 

a = 4.2 zb 0.1,q 0 = 0.2 zb 0.1 
qinf = 0.8 zb 0.03, r b = 6 zb 1 kpc 

4 

0 (-12758) 

0(25546) 

Incl. all substructures 

[Fe/H]-model 

f 0 = 0.6 zb 0.01, 7 = -0.15 zb 0.04 

2 



+ 





SPL 

a = 3.4 zb 0.1, q = 0.74 zb 0.02 

2 

-62 

110 

Einasto 

n = 2.3 zb 0.2, r e ff = 18 zb 1 kpc 
q = 0.77 zb 0.02 

3 

-24 

41 

BPL 

a in = 2.8 zb 0.1, a ou t — 4.3 zb 0.1 
r break = 29 zb 2kpc,g = 0.77 zb 0.02 

4 

-25 

52 

TPL 

ai = 2.8 zb 0.1, a 2 = 4.3 zb 0.2 
a 3 = 4.2 zb 0.5, q = 0.77 zb 0.02 
r breakl = 30kpc, r break2 = 55kpc 

4 

-26 

53 

[Fe/H]-model 

f 0 = 0.62 zb 0.01, 7 = -0.22 zb 0.06 

2 



+ 





Einasto-q(r) 

n = 6.1 zb 1.7, r e ff = 3 zb 2 kpc 
q b = 0.3 zb 0.05, qi n f = 0.9 zb 0.04 
r b = 9 zb 2 kpc 

5 

0 (-17675) 

6 

BPL-q(r) 

a in = 4.2 zb 0.3, a ou t = 4.5 zb 0.4 
r b reak = 32 zb 18 kpc ,qo = 0.3 zb 0.1 
qinf = 0.9 zb 0.04, r b = 9 zb 1 kpc 

6 

0 

12 

SPL-q(r) 

a = 4.4 zb 0.1,q 0 = 0.3 zb 0.1 
qinf = 0.9 zb 0.04, r b = 9 zb 1 kpc 

4 

-2 

0(35384) 


Notes. We give the type of model, the best-fitting parameters of the model, the number of free parameters, 
difference in log-likelihood from maximum likelihood value, and difference in BIC from minimum BIC value. 
Parameters that are kept fixed are highlighted in bold. 

a. A In Jz? max =ln Jz? max -max(ln Jz? max ), so AlnJz? max =0 means that the model has maximum likelihood. The 
value in parentheses is max(ln Jz? max ) • 

b. BIC — — 2 1nJz? ma;r + N p \n(Ndata)i an d A BIC=BIC — min {BIC), so the best model has AB IC=0. The 
value in parentheses is min(BIC). 



